*============================================================================================================================*
*** we conducted the following placebo tests and discuss the results in the paper: 
***   (1) we randomly assign the least‐cost distance values across all areas; 
***   (2) we randomly assign the piped water coverage values across all areas; 
***   (3) we randomly assign both the least‐cost distance and piped water coverage values across all areas.
*============================================================================================================================*


tempfile leastcostdist 

set more off 

*** set path here 

local regressors "l_pcincome lowwaterquality l_fisexp_pc l_fisexp_h_pc l_fisexp_inf_pc l_inst_inf_pc l_staff_inf_pc sharemanu rural2 mountain "   

*** Placebo tests 
** 1. randomly assign least cost distance 

use sample.dta, clear 

set obs 1000 
gen ivPlacebo=. 
gen ivtstat=. 
gen ivpvalue=. 
local k=1

forvalues i=1/1000 { 
	shufflevar l_cost_dist , cluster(prov_code) dropold	
	ren l_cost_dist_shuffled l_cost_dist  
	quietly xi: ivreg2 l_imr `regressors' (sharetapwater=l_cost_dist) i.prov_code    [w=totpop]  , partial(i.prov_code) robust 
	quietly replace ivPlacebo = _b[sharetapwater] in `k' 
	quietly replace ivtstat = _b[sharetapwater]/_se[sharetapwater] in `k' 
	quietly replace ivpvalue =2*ttail(e(df_m),abs(_b[sharetapwater]/_se[sharetapwater])) in `k' 
	
	local ++k 
} 

sum iv* 
winsor ivP, gen(ivPw) p(0.01)
duplicates drop ivPw, force  
hist ivPw, xtitle(Placebo IV Estimates) width(0.05) density  graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) saving(Placebo_IV_Estimates1, replace) 
graph export Placebo_IV_Estimates1winsor.jpg, replace  
hist ivPlacebo, xtitle(Placebo IV Estimates) width(0.05) density  graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) saving(Placebo_IV_Estimates1, replace) 
graph export Placebo_IV_Estimates1.jpg, replace  
hist ivt, xtitle(T-statistics of Placebo IV Estimates) width(0.1) density xline(-2 2) graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) saving(Tstat_Placebo_IV_Estimates1, replace)  
graph export Tstat_Placebo_IV_Estimates1.jpg, replace  
hist ivpv, xtitle(P-values of Placebo IV Estimates) width(0.025) density xline(0.05)  graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) saving(Pvalues_Placebo_IV_Estimates1, replace)  
graph export Pvalues_Placebo_IV_Estimates1.jpg, replace  

drop if ivPlacebo==. 
count if abs(ivt)>2 
count if abs(ivpv)<0.05 


** 2. randomly assign piped water coverage 
 
use sample.dta, clear 


set obs 1000 
gen ivPlacebo=. 
gen ivtstat=. 
gen ivpvalue=. 
local k=1

*** IV results 
forvalues i=1/1000 { 
	shufflevar sharetapwater , cluster(prov_code) dropold	
	ren sharetapwater_shuffled sharetapwater  
	quietly xi: ivreg2 l_imr `regressors' (sharetapwater=l_cost_dist) i.prov_code    [w=totpop]  , partial(i.prov_code) robust 
	quietly replace ivPlacebo = _b[sharetapwater] in `k' 
	quietly replace ivtstat = _b[sharetapwater]/_se[sharetapwater] in `k' 
	quietly replace ivpvalue =2*ttail(e(df_m),abs(_b[sharetapwater]/_se[sharetapwater])) in `k' 
	
	local ++k 
} 

sum iv* 
winsor ivP, gen(ivPw) p(0.01)
duplicates drop ivPw, force 
hist ivPw, xtitle(Placebo IV Estimates) width(0.05) density  graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) saving(Placebo_IV_Estimates2, replace) 
graph export Placebo_IV_Estimates2winsor.jpg, replace  
hist ivPlacebo, xtitle(Placebo IV Estimates) width(0.05) density  graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) saving(Placebo_IV_Estimates2, replace) 
graph export Placebo_IV_Estimates2.jpg, replace  
hist ivt, xtitle(T-statistics of Placebo IV Estimates) width(0.1) density xline(-2 2) graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) saving(Tstat_Placebo_IV_Estimates2, replace)  
graph export Tstat_Placebo_IV_Estimates2.jpg, replace  
hist ivpv, xtitle(P-values of Placebo IV Estimates) width(0.025) density xline(0.05)  graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) saving(Pvalues_Placebo_IV_Estimates2, replace)  
graph export Pvalues_Placebo_IV_Estimates2.jpg, replace  

drop if ivPlacebo==. 
count if abs(ivt)>2 
count if abs(ivpv)<0.05 


** 3. randomly assign both least cost distance and piped water coverage 

use sample.dta, clear 

set obs 1000 
gen ivPlacebo=. 
gen ivtstat=. 
gen ivpvalue=. 
local k=1

*** IV results 
forvalues i=1/1000 { 
	shufflevar sharetapwater l_cost_dist, cluster(prov_code) dropold	
	ren sharetapwater_shuffled sharetapwater  
	ren l_cost_dist_shuffled l_cost_dist  
	quietly xi: ivreg2 l_imr `regressors' (sharetapwater=l_cost_dist) i.prov_code    [w=totpop]  , partial(i.prov_code) robust 
	quietly replace ivPlacebo = _b[sharetapwater] in `k' 
	quietly replace ivtstat = _b[sharetapwater]/_se[sharetapwater] in `k' 
	quietly replace ivpvalue =2*ttail(e(df_m),abs(_b[sharetapwater]/_se[sharetapwater])) in `k' 
	
	local ++k 
} 

sum iv* 
winsor ivP, gen(ivPw) p(0.01) 
duplicates drop ivPw, force 
hist ivPw, xtitle(Placebo IV Estimates) width(0.05) density  graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) saving(Placebo_IV_Estimates3, replace) 
graph export Placebo_IV_Estimates3winsor.jpg, replace  
hist ivPlacebo, xtitle(Placebo IV Estimates) width(0.05) density  graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) saving(Placebo_IV_Estimates3, replace) 
graph export Placebo_IV_Estimates3.jpg, replace  
hist ivt, xtitle(T-statistics of Placebo IV Estimates) width(0.1) density xline(-2 2) graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) saving(Tstat_Placebo_IV_Estimates3, replace)  
graph export Tstat_Placebo_IV_Estimates3.jpg, replace  
hist ivpv, xtitle(P-values of Placebo IV Estimates) width(0.025) density xline(0.05)  graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) saving(Pvalues_Placebo_IV_Estimates3, replace)  
graph export Pvalues_Placebo_IV_Estimates3.jpg, replace  

drop if ivPlacebo==. 
count if abs(ivt)>2 
count if abs(ivpv)<0.05 



** 4. randomly assigned the water pollution level and the corresponding least-cost distance to IMR locations. Then, based on the randomly assigned water pollution levels, we divide the sample by the pollution level and re-estimate the full model (Column 7 in Table 4). 

use sample.dta, clear  

gen prov_highpol=prov_code   
recode prov_highpol (43=42) (64=63)  

set obs 1000 
gen dif=. 
gen pvalue=. 
gen ivplacebo1=. 
gen ivplacebo2=. 

local k=1

*** IV results 
forvalues i=1/1000 { 

	shufflevar waterquality l_cost_dist, cluster(prov_code) dropold	
	ren waterquality_shuffled waterquality  
	ren l_cost_dist_shuffled l_cost_dist  
	
	xtile waterquality_group2=waterquality, n(2) 

	quietly xi: ivreg2 l_imr `regressors' (sharetapwater=l_cost_dist) i.prov_code    [w=totpop]   if waterquality_group2==1, robust partial(i.prov_code) 
	quietly replace ivplacebo1 = _b[sharetapwater] in `k' 

	quietly xi: ivreg2 l_imr `regressors' (sharetapwater=l_cost_dist) i.prov_code    [w=totpop]   if waterquality_group2==2, robust partial(i.prov_code) 
	quietly replace ivplacebo2 = _b[sharetapwater] in `k' 


	gen cons=1
	quiet ivreg2 l_imr c.(l_pcincome lowwaterquality l_fisexp_pc l_fisexp_h_pc l_fisexp_inf_pc l_inst_inf_pc l_staff_inf_pc sharemanu rural2 mountain )#i.waterquality_group2 i.prov_code#i.waterquality_group2 c.cons#i.waterquality_group2 (c.sharetapwater#i.waterquality_group2 = c.l_cost_dist#i.waterquality_group2), nocons robust partial(i.prov_code#i.waterquality_group2)   
	
	quietly replace dif = 1.waterquality_group2#c.sharetapwater - 2.waterquality_group2#c.sharetapwater in `k' 
	quietly test 1.waterquality_group2#c.sharetapwater = 2.waterquality_group2#c.sharetapwater
	quietly replace pvalue = r(p) in `k' 

	drop waterquality_group2 cons 
	local ++k 
} 

sum dif pvalue 

winsor dif, gen(difw) p(0.01) 
* duplicates drop difw, force 
winsor ivplacebo1, gen(ivplacebo1w) p(0.01) 
winsor ivplacebo2, gen(ivplacebo2w) p(0.01) 
duplicates drop ivplacebo1w, force 
duplicates drop ivplacebo2w, force 
hist ivplacebo1w, graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255))  xtitle(Placebo IV Estimates (Low Pollution Areas)) width(0.05) density  
graph export Placebo_IV_Estimates_LowP.jpg, replace  
hist ivplacebo2w, graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255))  xtitle(Placebo IV Estimates (High Pollution Areas)) width(0.05) density  
graph export Placebo_IV_Estimates_HighP.jpg, replace  

twoway (histogram ivplacebo1w, width(0.2) density)   (histogram ivplacebo2w, width(0.2) density fcolor(none) lcolor(black)), graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255)) legend(order(1 "Low Pollution Areas" 2 "High Pollution Areas" )) xtitle(Placebo IV Estimates) 
graph export Placebo_IV_Estimates_2hist_pollution.jpg, replace  

hist difw, graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255))  xtitle(Differences in IV Estimates of Low- and High-Pollution Areas) width(0.05) density 
graph export Placebo_LowHighPDif_Estimates_winsor.jpg, replace  
hist dif, graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255))  xtitle(Differences in IV Estimates of Low- and High-Pollution Areas) width(0.05) density xline(-0.311) 
graph export Placebo_LowHighPDif_Estimates.jpg, replace  
hist pvalue, graphregion(fcolor(255 255 255)) graphregion(lcolor(255 255 255))  xtitle(P-values of Differences in IV Estimates of Low- and High-Pollution Areas) width(0.05) density xline(0.05)  
graph export Placebo_LowHighPDif_pvalues.jpg, replace  

